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Abstract 

A new code for the scale evolution of modified-minimal-subtraction-scheme 
parton densities is described. Through next-to-leading order the program uses 
exact splitting functions. In next-to-next-to-leading order approximate split- 
ting functions are used. For efficiency the program includes analytical results 
for the evaluation of the weights required for the integrations over the lon- 
gitudinal momentum fractions of the partons. It also incorporates the opera- 
tor matrix elements required for the matching conditions across heavy flavor 
thresholds in higher order perturbation theory. The more efficient handling of 
the weights implies that the code is faster than similar evolution codes in all 
modes of operation. The program is written in the C programming language. 

PACS numbers: ll.lOJj, 12.38Bx, 13.60Hb, 13.87Ce. 
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1 Program Summary 



Title of Program: ADENS 

Computer: AlphaStation 4/500 

Operating system: OSF V3.2 

Programming language: C 

Number of lines in distributed program: 11658 

Keywords: parton density, evolution, numerical solution, splitting function, 
next-to-next-to-leading order 

Nature of physical problem: solution of the parton density evolution equations 
with LO, NLO and NNLO splitting functions and NLO, NNLO heavy flavor 
threshold matching conditions 

Method of solution: x-space integration with analytic evaluation of weights 
Typical running time: see table in Section 7 
The program www site: 

http: / /insti. physics. sunysb.edu/ ~chuvakin/ adens-24.0.tar.gz 
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http://msti.physics.sunysb.edU/~smith/adens-24.0.tar.gz 
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Contents 



2 Introduction 



Deep-inelastic lepton-liadron scattering experiments probe the internal struc- 
ture of hadrons. The lepton-hadron inclusive cross sections may be written 
in terms of structure functions, which depend on the virtuality of the probe 
Q"^ . Three structure functions Fi, and are necessary to describe neu- 
tral current (photon and Z-boson exchange) and charged current (VT-boson 
exchange) reactions. In perturbative quantum chromodynamics (pQCD) the 
probe interacts with partonic constituents of the hadron. There are probability 
densities /(x,/i^) to find partons carrying a fraction a; (0 < a; < 1) of the lon- 
gitudinal momentum of the hadron at a mass factorization scale yU. Therefore 
the Fj, i = 1,2, L also depend on x and fi. 

The operator product expansion (OPE) allows the structure functions to be 
written as convolutions of the parton (quark and gluon) probability densities 
with partonic hard scattering cross sections (or coefficient functions). The lat- 
ter can be calculated in pQCD. Even though the former cannot be calculated in 
pQCD, their /i dependence is determined by a set of integro-differential equa- 
tions, the (Dokshitzer-Gribov-Lipatov)-Altarelli-Parisi evolution equations 
which follow from renormalization group analysis. Discussions of the pQCD 
description of deep inelastic scattering reactions are available in and 0. 
The probability densities and splitting functions are defined in the modified- 
minimal-subtraction (MS) scheme. 



For simplicity we will call the above equations the evolution equations. They 
describe processes where a massless parton (quark or gluon) carrying a fraction 
of the longitudinal momentum of the incoming hadron radiates a massless 
parton and becomes a (different) massless parton with a different momentum 
fraction. The probability for this process to happen is determined by splitting 
functions which are computed order-by-order in pQCD. The leading-order 
(LO) and next-to-leading order (NLO) splitting functions have been known 
for some time [^, [^], [0, [0, [0, [0 and the results are summarized in a 
convenient form in |^ . Recently some moments of the next-to-next-to-leading 
order (NNLO) splitting functions have been calculated in [|lOl, see also |ll] 



and . If the x-dependence of the quark and gluon densities in a hadron are 
parametrized at one value of fi, (say at /io,) then the solutions of the evolution 
equations with the above LO, NLO or NNLO splitting functions yield the 
X dependence of the massless parton densities at a different fi. There is a 
second scale in the pQCD theory, the renormalization scale, which appears in 
argument of the the running coupling a^. It is usually set to be the same as 
the mass factorization scale fi so ag = as(/x^). 
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The flavor dependence of the quark and anti-quark densities is governed by 
the flavor group, which is SU(2) for the up and down quarks, SU(3) for the 
up, down and strange quarks etc. Therefore is is convenient to form flavor 
non-singlet and flavor (pure) singlet combinations of densities. The former 
have their own evolution equations. The latter mix with the gluons and the 
combined evolution is described by matrices which obey coupled integro- 
differential equations. 



A number of methods to solve the evolution equations for the parton densities 
have been proposed, including direct x-space methods, |13|, [|T^, [|T5|, ||T6[, |T7[] , 
TSj , orthogonal polynomial methods |]19|, [^, and Mellin-transform methods 
23]. A compilation of parton density sets is available in p3[ . 



21 



The best method, which should be both accurate and fast, depends on region 
chosen in x and /x^. Currently the requirements are that the code be able to 
evolve densities from a minimum /x^ near 0.26 GeV^ up to a maximum fx^ near 
10^ GeV^ required for QCD studies for the future Large Hadron Collider at 
CERN. The range in x is from a minimum value near 10~^ up to a maximum 
near unity. We use the direct x-space method, with the following additional 
features. 



One of our aims is a better treatment of parton density evolution for "light" 
u, d and s quarks near the heavy flavor thresholds chosen to be at the charm 
and the bottom quark masses {m^. and respectively). The parton density 
description must be modifled to incorporate new c and h "heavy" quark den- 
sities as the evolution scale increases. The implementation of the NLO and 
NNLO matching conditions across heavy flavor thresholds in the variable fla- 
vor number schemes (VFNS) [^, |2^, [^, ||2^ involve large cancellations 



between various terms in the expressions for the structure functions. Poor nu- 
merical accuracy in the solution for the evolution of the parton densities at 
small scales would spoil these cancellations and ruin the VFNS predictions. 
We achieve the required accuracy by avoiding one numerical integration in 
our program so we analytically calculate the weights for the exact LO, the 
exact NLO and the approximate NNLO splitting functions. The approximate 
NNLO splitting functions are taken from p8[,|2^, while the relevant opera- 



tor matrix elements (OMEs), which provide the matching conditions on the 



parton densities across heavy flavor thresholds, are taken from |30 



Since we start the scale evolution from a set of densities (input boundary 
conditions) at a low scale = /xq <^ w^c, the running coupling ^^(/i^) is large. 
We therefore use the exact solution of the NLO equation for and match the 
values on both sides of the heavy flavor thresholds to three decimal places. We 
mention here that the NNLO matching conditions on as across heavy flavor 
thresholds are available in ||31| and Our program evolves both light and 



heavy parton densities in LO, NLO and NNLO from a minimum x equal to 
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10~^ to a maximum x equal to unity, a mimimum = 0.26 (GeV)^ in LO 
and /i^ = 0.40 (GeV)^ in NLO and NNLO and and a maximum /i^ = lO'^ 
(GeV)^. Results have been published in 
a detailed write up of the program. 



26| , [p7[| and |^3[. Here we give 
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3 The evolution equations 



3.1 Definitions of densities 



We evolve combinations of up (m). down (rf), strange (s), charm (c) and bot- 
tom (6) quark densities whicli transform appropriately under the flavor group. 
Hence we define flavor-non-singlet valence quark densities by 

fk-kirif, X, //^) = fkirif, X, iJ?) - fk{nf, x, /j,'^) , k ^ u,d . (3.1) 



The flavor-singlet quark densities 

fq{nf: X, fJ,^) = fk+ki'^f: X, //^) (3.2) 

k=l 



are deflned in terms of the expression 

fk+k{rif, X, ji^) = fk{nf, X, 11^) + /fc(n/, x, i^"^) , k = u, d, s,c,b, (3.3) 

when Uf — 5. Then the flavor-non-singlet sea quark densities are 

fg^iuf, X, 11^) = fk+-k{nf, X, li^) - —fqirif, X, //^) . (3.4) 

These equations will be discussed further in the next section. 

3.2 The evolution equations 



A typical evolution equation is that for a flavor-non-singlet parton density 

r^(x,/.^) 

9 rNSr 2^ (^s{lJ^) } dx y 



din fj? 



rNS/ 2\ ^s[fJ' ) f ^X ,y 2\ rNS/ 2\ /o 



y 



where P^^{y/x, ii'^) is a non-singlet splitting function, and q;s(/x^) is the run- 
ning coupling. 
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The sphtting functions in the evolution equations can be expanded in a per- 
turbation series in into LO, NLO and NNLO terms as follows 

P{z,i,') = P^'\z,i,') + (^)PW(.,/.^) + (^)^P(^)(.,/.^). (3.6) 



The non-singlet combinations of the qriq,) to qs{qs) splitting functions, where 
the subscripts r, s denote the flavors of the (anti)quarks q and q respectively 
and satisfy r,s = l,---,ny, can be further decomposed into flavor diagonal 
parts proportional to Srs and flavor independent parts. In LO there is only 
one non-singlet splitting function Pgg but in NLO it is convenient to form two 
combinations from P„„ and P„^ as follows 



P+ Pqq + Pqq i 

P- = Pm-Pm- (3-7) 

These splitting functions are used to evolve two independent types of non- 
singlet densities, which will be called plus and minus respectively. They are 
given by 

^ fk--k{nf,x,fJ,^) . (3.8) 

Since the general formulae in Eqs. (3.1)- (3.4) are rather involved the easiest 
way to explain the indices is by explicitly giving the combinations we use. For 
j — 1,2 we have 

f-^u-u,f2^d-d, (3.9) 



which are used for all flavor density sets. Then for three-flavor densities i — 
1,2,3 and we deflne 

ft^u + u- E(3)/3 , f+^d + d- E(3)/3 , 

/3+ = s + s-E(3)/3, (3.10) 

where E(3) = fq{S) — u + u + d + d + s + s. These densities should be used 
for scales fj, < rric. For four-flavor densities i = 1, 2, 3, 4 and we deflne 

f+^u + u- E(4)/4 , f+^d + d- E(4)/4 , 

f+^s + s- E(4)/4 , f+^c + c- E(4)/4 , (3.11) 

where E(4) = /^(4) = c + c + E(3). These should be used for scales in the 
region rric < < rrib- For flve-flavor densites i = 1, 2, 3, 4, 5 and we deflne 
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/+ = ii + iZ-E(5)/5, 

/3+ = ,9 + S-S(5)/5, 



/2+ = d + J-S(5)/5, 
/4+ = c + c-E(5)/5, 



(3.12) 



where E(5) = fq{5) — b + b + S(4). These should be used for scales fx > nib. 

If we define t = ln(//7(l GeV^) then we need to solve the four evolution 
equations 



dft{y,t) 


a sit) 


dt 


27r 


dfiiy,t) 


_ (^s{t) 


dt 




dfs{y,t) 


_ <^s{t) 


dt 


27r 




_ Ois{t) 


dt 


2tt 



-P4-,t)ft{x,t), 

X X 



-P-{'-,t)frix,t), 
J X X ■' 



I - \p<A-.t)flM + P99{-.t)f!M 



(3.13) 



(3.14) 



(3.15) 



(3.16) 



where for ji < rric wc set i = 1,2,3, j = 1,2, = S(3) and the gluon is 



a three-flavor gluon. When rric < l^i < rrib, we use i = 1,2,3,4, j = 1,2, 



= E(4) and the gluon is a four-flavor gluon. Finally when > mf,, we set 



1 — 1,2, 3, 4, 5, j = 1, 2, = E(5) and the gluon is a five-flavor gluon. Note 
that since NNLO splitting functions are approximate we provide the high and 
low estimate for each splitting functions labeled A and B. For all calculations 
we use their average so that the error is minimized. 

The densities satisfy the momentum conservation sum rule which we write in 
terms of the u, d, ..b (anti)-quark and gluon densities as 



dxx 



u{x, fx^) + u{x, 11^) + d{x, 11^) + d{x, /x^ 



+s{x, jj?) + s{x, /i^) + [c{x, jj?) + c{x, iJ^)]9{iJi'^ - ml) 
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+ [b{x, yU^) + b{x, ^^)\9{^^ - ml) + g{x, 



1. 



(3.17) 



Also the quark constituents carry all the charge, isospin, strangeness, charm 
and bottom quantum numbers of the nucleon so they also satisfy the other 
standard sum rules for the conservation of these quantities, see @, 0. 
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4 Direct x-space method of solution and initial conditions 



4.1 The method 



Our choice of the direct x-space method is motivated by the necessity to step 
densities across heavy flavor thresholds using LO, NLO and NNLO boundary 
conditions. The procedure of doing this with Mellin moments would involve 
taking moments of the densities and then inverting moments several times. 
The direct x-space method is much more intuitive and straightforward. The 
main features of this method are linear interpolation over a grid in x and 
second-order interpolation over a grid in t. Let us describe it in more detail 
to point out where we differ from the method in 



First we consider the x-variable in the evolution and write the right-hand-side 
of the evolution equation Eg. (^.5|) for the non-singlet density as 

JM=/^£^£pf£5),W, (4.1) 



where Xq < a; < 1 , 

q{x) = xf{x) , (4.2) 

and 

Xo<Xi< ... <Xn< Xn+l = 1 , (4.3) 

with q{xn+i) = q{l) = 0. Between grid points Xj and Xj+i, x is chosen so that 

q{x) = (1 - x)q{xi) + xq{xi+i) , (4.4) 

with X = {x — Xi)/{xi+i — Xi). Using this relation we convert the integral into 
a sum 

n+l 

K^o) = ^(^»' ^o)qix,i) , (4.5) 

where the weights are (in all orders LO, NLO and NNLO) 

w{xo,xo) = Si{si,so) 

w{xi, xo) = 5*1(5^+1, Si) - S2isi, Si_i) , (4.6) 
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where Sj = x^/xi and 



S,{u,v) = ^ l{z-u)P{z)- 



V — u J ' ' ■ ■ z ' 



lit d 7 

S2{u, v) = / {z - v)P{z)— . (4.7) 

V — u J z 



u 



In the above formula P{z) denotes the sphtting function of the correspond- 
ing order in and type (non-singlet, singlet, etc.) We use the LO and NLO 
splitting functions in and the approximations to the NNLO splitting func- 



tions from [^] and For completeness the latter are given in Appendix A. 
We have calculated the integrals in Eq. ([4.7|) analytically and the results are 
in the computer program. This yields the formula in Eq. (|4.5| ) describing the 
grid for the x variable. Note that the weights w^'^\ w^^^ and w^'^^ are those for 
the exact LO, the exact NLO and the approximate NNLO splitting functions 
respectively. Thus, for the singlet case, we have 



d{xoT.{xo)) (0)/ N , "s fi)/ \ , /"s 



X XjS(xj) 

+ [wfJ{Xi,Xo) + ^wlj]^{Xi,Xo) + {^fw^qJ{Xi,Xo)] 

X Xig{xi) , (4. 



where S is either S(3), S(4) or S(5) depending on the scale. 

Now consider the variation in the variable t. For each Xi we pick a grid in 
t labelled by distinct points tj. Then, for example, the non-singlet equation 
becomes 



q'{xi,tj) = ^^4~^Y1 [w'i\xk,Xi) + ^^^^w'i\xk,Xi) 
271 27r 

+ (^)2^(?)^^^^^^^^^^^^^^^ (4.9) 

where q'{xi,tj) denotes the derivative with respect to t evaluated at t = tj. In 
compact notation this equation can be rewritten as 

q'j = wqj + S, (4.10) 
with 5* being the sum of the terms on the right hand side of Eq. (|4.9|) excluding 
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the j-th term. 



For t between the grid points and tj we interpolate the parton density 
using quadratic interpolation as follows: 

q{xi,t) = at'^ + bt + c. (4.11) 



Thus we relate the value of q at the point tj to that of q at the point by 

q{xi,tj) = q{xi,tj_i) + ^[q'{xi,tj) + q {xi,tj-i)]Atj , (4.12) 



where Atj = tj — tj-i- This equation can also be written more compactly as 

qj = q,-i + lii-i + qj)^tj. (4.13) 



The resulting system of two linear equations in Eq. (^4.10|) and Eq. ( [4.13|) for 
qj and has the solution 

2a-i + (^-. + S)At, 

2-wAtj ^ ' 



Then we find q^ from Eq. ( [4.10D . Applying the same procedure to the gluon 
and singlet combinations involves four equations because we have to compute 
both the densities and their derivatives. 

The evolution proceeds from the initial /Xq = /^lo (^^ A^o — /^nlo) ^'^ ^^^t 
heavy flavor threshold at the scale /i^ = m^. Next the charm density is intro- 
duced in NNLO (a^-order terms) and all the four-flavor densities are evolved 
from the new boundary conditions in Section 4.2. This evolution continues 
up to the transition point /x^ = m^, where the same procedure is applied to 
generate the bottom quark density. From that matching point all five-flavor 
densities are evolved up to all higher /i^ scales starting from the boundary 
conditions in Appendix B. 

Since the weights for the calculation are computed analytically from the LO, 
NLO [0 and NNLO (|2|],[gg) MS splitting functions we remove possible 



instabilities in the numerical integrations. Hence the program is very efficient 
and fast. The results from the evolution code have been thoroughly checked 
against the tables in the HERA report [jl6[ and they agree to all five decimal 
places. 
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4-2 The initial conditions 

The GRV98 three-flavor LO and NLO parton density sets contain in- 
put formulae at low scales fi < rUc which are ideal as initial values for our 
parametrizations. Therefore we start our LO evolution using the following 
input at /ig = ij,Iq = 0.26 GeV^ 



= 1.239 x^-^^ (1 - xf^^ (1 - l.SVi + 9.5a;) 

= 0.614 (l-x)°-^ xu^{x,^ilo) 
a;(/j(3, X, /io) - /s(3, X, /Xq)) = xA{x, /Xlo) 

= 0.23 a;°-^^ (1 - a;)"-^ (1 - 12.0Vx + 50.9a;) 
a;(/d(3, X, /io) + /ii(3, X, /Xq)) = x{u + d){x, /iLo) 

= 1.52 a;°-i^ (1 - xf-^ (1 - 3.6^^ + 7.8x) 
xfg{3,x,nl) =xg{x,filo) 

= 17A7 x^-^{l-xf-^ 

^/s(3, X, /Xq) = x/g(3, X, /ig) = Xs{x, /iLo) 

= xs(x,/i2o) = 0. (4.15) 

Here A = d — u is used to construct the non-singlet combination. 

We start the corresponding NLO evolution using the following GRV98 input 
at fil = /i^LO = 0-40 GeV^ 



X X, /ig) — XU-ui^X, /i^Lo) 

= 0.632 x" ''^ (1 - a;)=' °^ (1 + 18.2a;) 

^fd-di'^j ^) f'o) — xd^(x, /iNLo) 

= 0.624 (l-x)^-° XM,(x,/i^Lo) 
^(/d(3, a;, /ig) — /u(3, x, /ig)) = xA(x, /i^Lo) 

= 0.20 a;°-^^ (1 - x)^^'"^ (1 - 13.3Vx + 60.0a;) 

a;(/<j(3, X, /ig) + /s(3, X, /ig)) = X(U + J)(X, /iNLo) 

= 1.24 x°-^° (1 - x)*^-^ (1 - 2.3v^ + 5.7x) 
^/g(3, X, /ig) = xg(x, /i^Lo) 

= 20.80 x^-^(l-x)^-^ 

a^/s(3, X, /ig) = x/g(3, X, /ig) = Xs(x, /iNLo) 

= xs(x,/i^Lo) = 0. (4.16) 
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We start the corresponding NNLO evolution using the same NLO input and 
starting scale as above. 



4-3 The calculation of the running coupling 



The heavy quark masses nic = 1.4 GeV^, ruh = 4.5 GeV^ are used throughout 
the calculation. We also use the exact solution (as opposed to a perturbative 
solution in inverse powers of ln(/i^/A^)) of the differential equation 



d ln(;u2) 



IGvr^ 



(4.17) 



for the running coupling ^^(/i^). Here (3q = 11 — 2nj/3 and Pi = 102 — 38r2y/3. 
Another way of writing this equation is 



In 



V*-EXACT 



f> y 



PI 



In 



An 



+ 



(4.18) 



The values for A^xact carefully chosen to obtain accurate matching of 
at the scales ml and m^. We used the values Aj^xact = 299.4, 246, 167.7, 67.8 
MeV/c^ respectively in the exact formula (which yields a 



EXACT/ 



EXACT ^ 

(3,4,5,6) 
LO 



mt 



0.205, a 



EXACT ^ 



0.319, «f^ACT(^ 



NLO J 



= 0.114, 
0^578 ) and 



204, 175, 132,66.5 MeV/c^ respectively (which yields a 



LO, 



0.125, a^^iml) 



0.232. 



0.362, a^O(^^ 



Lo) 



0.763 ) for the 



LO formula (where /3i = 0). There is a NNLO discontinuity of aproximately 
two parts in one thousand in the running coupling across heavy flavor thresh- 
olds We have ignored this effect to focus on the numerically more 



significant matching of the flavor densities. 



4.4 The evolution process 



Three flavor evolution proceeds from the initial /^g to the scale /x^ 
(GeV^)2. At this point the charm density is then defined by 



1.96 



/e+c(% + 1, ml) = alirif, ml) [AHil) ® f^{nf, ml) 
+4.(l)®/'(^/'^c)], 



(4.19) 



with Uf = 3 and = as/ An. We have suppressed the x dependence to make 
the notation more compact. The symbol denotes the convolution integral 
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/®^? = I fix/y)g{y)dy/y, where x < y < 1. The OME's ^^^(/iVm^), 
Agglfi"^ /m^) are given in |0 and are also listed in Appendix B. The rea- 
son for choosing the matching scale fi at the mass of the charm quark rric is 
that all the ln(/i^/m^) terms in the OME's vanish at this point leaving only 
the nonlogarithmic pieces in the order OME's to contribute to the right- 
hand-side of Eq.(4.19). Hence the LO and NLO charm densities vanish at the 
scale fi = rric- The NNLO charm density starts off with a finite x-dependent 
shape in order a^. Note that we then order the terms on the right-hand-side 
of Eq. ( |4.19| ) so that it contains a product of NLO OME's and LO parton 
densities. The result is then of order and should be multiplied by order a° 
coefficient functions when forming the deep inelastic structure functions. 

The four-flavor gluon density is also generated at the matching point in the 
same way. At fj, = rric we define 



+al{nf,ml] 



(4.20) 



The OME's A^g^qifJ''^ /ml), A^g^qifJ^'^/ml) are given in |^ and are also fisted 
in the Appendix B. The four-flavor light quark (u,d,s) densities are generated 
using 



,2\ 



+af(n;,m,^)A^;^Q(l) ® • (4-21) 

The OME ^^^^^(/iVm^) is given in (as well as in Appendix B) and the 
total four-flavor singlet quark density folows from the sum of Eqs. ( |4.19|) and 
( irap . In Eqs. ( |CT| ) and (^^) we set n/ = 3. The remarks after Eq. ( ^:TP| ) 



are relevant here too. 



Next the resulting four-flavor densities are evolved using the four-flavor weights 
in either LO, NLO and NNLO up to the scale = ml = 20.25 (GeV^)^. The 
bottom quark density is then generated at this point using 



|(S) 



(4.22) 



and the gluon and light quark densities (which now include charm) are gener- 
ated using Eqs.( [4.19D -( [4.21D with n/ = 4 and replacing by m^. Therefore 
only the nonlogarithmic terms in the order OME's contribute to the match- 
ing conditions on the bottom quark density. Then all the densities are evolved 
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up to higher five-flavor set with either LO, NLO and NNLO sphtting 

functions. This is vahd until /i = ^ 175 GeV^ above which one should 
switch to a six-flavor set. We do not implement this step because the top 
quark density would be extremely small. 

The procedure outlined above generates a full set of parton densities (gluon, 
singlet, non-singlet light and heavy quark densities,) for any x and yU^ from the 
three-flavor LO, NLO and NNLO inputs in Eqs.( P^ and pi^ ). Note that 



we have used /i^ for the factorization and renormalization scales in the above 
discussion. In the computer program we use the notation that denotes 
these scales, since this is done in all the previous computer codes for the 
parton densities. 
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5 Input peirameter description and usage 



To prepare the program for use unpack the distribution package adens-24.tar.gz 
by typing tar -xzf adens-24.tar.gz. The resulting directory will contain the fol- 
lowing files 



head.h 
main.h 
main . c 
1-a-w . c 
nl-a-w. c 
alpha . c 
init . c 
polylo . c 
intpol . c 
evolver . c 
thresh. c 
a-coef s . c 
loader . c 
quadrat . c 
daind . c 
integrands . c 
grids . c 
weights . c 
nnl-a-w. c 
wgplg.c 

evolution_parameters . input 

makefile 
my_howto .tex 
sample . out 



To build the executable on a machine with a gcc compiler type make . The 
executable named adens.x will be produced. To run the code just run the file 
adens.x. Some debugging information may appear on the standard output. 



Here is the parameter file (evolution_parameters. input) explanation with 
default values shown: 
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0.204e0 


LambdaLO-3 


LO A for A^. =3 


0.175e0 


LambdaLO-4 


LO A for Nf =4 


0.132e0 


LambdaLO-5 


LO A for Nf =5 


0.306e0 


LambdaNL03 


NLO A for Nf =3 


0.257e0 


LambdaNL04 


NLO A for Nf =4 


0.1734e0 


LambdaNLOS 


NLO A for Nf =5 


0.2994e0 


LambdaENL03 


Exact A for Nf =3 


0.246e0 


LambdaENL04 


Exact A for A^^ =3 


0.1677e0 


LambdaENLOS 


Exact A for Nf =3 


0.40e0 


Qinitial2 


Initial to start evolution 


1.96c0 


QcharmMass 


Mass of first heavy quark c 


20.25e0 


QbottomMass 


Mass of second heavy quark b 


1.96e0 


QcharmThreshold 


Charm threshold 


1.96e0 


AlphaCharmThreshold 


C threshold used for ag 


20.25c0 


QbottomThreshold 


Bottom threshold 


20.25e0 


AlphaQbottomThreshold 


B threshold used for as 


lOOO.OeO 


Qfinal2 


Final g2 


130 


tGridSize 


grid size 


200 


xGridSize 


X grid size 


130 


xGridSplit 


X split between log and linear 


l.Oe-5 


xinitial 


X initial 


0.2e0 


xS])lit 


,r at tlie s])lit htw log aiid linear 


l.OOeO 


xFinal 


X final (always 1) 





DebugLevel 


Error message detail (0-5) 


1 


GraphVsX 


Plotting data files are versus x {1) 

or (0) 


1 


Order 


LO/NLO/NNLO for 0,1,2 





DoFortran 


Produce (1) or no (0) data files for 
CSN/BMSN Fortran programs 
(1-yes, 0-no) 


1 


AlphaDoSeparateThreshold 


Use separate thresholds for ag (1- 
yes, 0-no) 
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1 


AlphaUseExact 


Use exact GRV98-style (1-yes, 
0-no) 





ThreeFlavorMode 


Calculate GRV98-style densities 
with no heavy flavors (1-yes, 0- 
no) 





GraphAU 


riot all data points (1-yes, U-no) 





NNLOmultiOrderCHARM 


Use our proper order NNLO 
heavy flavors (1-yes, 0-no) 


1 


DoBottomThreshold 


Generate bottom (1-yes, 0-no) 





LoadWeightsMadeBefore 


Use ready weights if available (1- 
yes, 0-no) 


1 


DoNotDump Weights 


Dump weight for future use as the 
option above (1-yes, 0-no) 





NL04NNL0 


Use NLO weights for NNLO cal- 
culation (1-yes, 0-no) 



The flrst set of Lambdas are used for LO calculations. The second set are used 
for NLO and NNLO calculations if the exact as is not requested (AlphaUse- 
Exact=0). The next set (LambdaENL03, LambdaENL04, LambdaENLOS ) 
are used for the exact solution of the differential equation for as as proposed 



in the GRV98 paper |2^. The code that calculates the exact might use its 
own set of flavor thresholds (which means that the number of flavors used for 
as can be reset independently from the regular heavy flavor threshold as done 
in m). 



Next we give the limits and the heavy masses: the initial and flnal Q^, the 
charm and bottom masses (used in threshold calculations), the heavy flavor 
thresholds and the separate as thresholds. Next follow the grid sizes in x and 
together with x initial and flnal (always 1) and the switch point between 
logarithmic and linear grids in the x dimension. The x grid always starts as 
logarithmic and then becomes linear at higher x, usually at a value of the 
order of 0.1 (xGridSplit parameter). 



The last group of parameters contains various control values that set the modes 
of the computation: 

DebugLevel , controls the amount of generated error, warning and information 
messages. 

Graph VsX , controls the printing of the output data for plotting (flrst column 
is either x or Q^, then subsequent columns will contain density values for var- 
ious or x), 
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Order, sets calculation order (use 0,1,2 for LO,NLO,NNLO), 

DoFortran, sets whether to dump interpolated densities on a special grid for 

future use in Fortran code for the calculation of structure functions ; CSN and 



BMSN refer to VFNS schemes which are explained in |p5|| , 
AlphaDoSeparateThreshold, sets whether we use a separate threshold for as 
(used, for instance for GRV98 set where Uf for densities is always 3 and Uf 
for as goes from 3 to 5, 

AlphallseExact, sets whether to use exact (differential equation solution) as 
for NLO and NNLO calculation, 

ThreeFlavorMode, sets whether to run GRV98 mode (no heavy flavors, Uf = 3 
for all Q^), 

GraphAll, controls the amount of graphing and printing output (either all data 
points or the special grid defined in the file main.h, that contains some fa- 
vorite values (for more see Section 7)), 

NNLOmultiOrderCHARM, activates NNLO threshold calculation using proper 
order combinations (this mode requires one to first run the LO and NLO cal- 
culations), 

DoBottomThreshold, enables the bottom density, 

LoadWeightsMadeBefore, turns on and off the loading of weights computed in 
the prior runs, 

DoNotDump Weights, sets whether to save computed weights to disk for future 
use, 

NLO4NNLO, sets whether NLO weights are used for the NNLO calculation 
(thus having only the boundary condition in NNLO). 

Some common parameter settings and typical grid sizes for popular evolutions 
are shown in Section 7. 
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6 Description of the program 



6.1 Program module summary 



main.c 


The main program, input and output 


1-a-w.c 


Calculation of LO weights 


nl-a-w.c 


Calculation of NLO weights 


alpha, c 


Calculation of as 


init.c 


Definition of initial functions 


polylo.c 


Calculation of polylogarithms 


intpol.c 


Interpolation routine 


evolver.c 


Evolution process subroutine 


thresh. c 


Threshold handling subroutine 


a-coefs.c 


OMEs for thresholds 


loader.c 


Datafile reading subroutine 


quadrat. c 


Gaussian integration subroutine 


daind.c 


Another integration subroutine 


integrands, c 


Heavy fiavor integrand calculation routine 


grids. c 


Grid generation routine and 




memory management routines 


weights.c 


Weight table handhng routine 


nnl-a-w.c 


Calculation of NNLO weights 


wgplg.c 


Calculation of high order polylogarithms 



6.2 main.c 

subroutines: 
none. 

The main program module contains input handling from the parameter file, pa- 
rameter verification, calls to grid generating routines (MakeXGrid, McikeT- 
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Grid), resets for all density arrays (array q) and their derivatives (array qp). 
It also includes calls to the generation of weights (analowgts, ananlowgts 
), the calls to evolution and threshold routines (evolver, threshold) that do 
the actual work. Also it contains some pre-output density processing and the 
results provided in various formats for both viewing and plotting. 

6.3 l-a-w.c 

subroutines: 

int analowgts(int nf,int loadWgts), 
int computeLOwgts(int nf), 
double sqq(double x,double y), 
double sgg (double x,double y). 

Analytically computes or reads from the file the LO weights for the evolution 
equations. 

6.4 nl-a-w.c 

subroutines: 

int ananlowgts(int nf,int loadWgts), 
int computeNLOwgts(int nf), 
double slff (double x, double y, int nf), 
double s2ff (double x, double y, int nf), 
double slfg(double x,double y, int nf), 
double s2fg(double x,double y, int nf), 
double slgf (double x,double y, int nf), 
double s2gf(double x,double y, int nf), 
double slgg(double x, double y, int nf), 
double s2gg(double x, double y, int nf), 
double slff_plus(double x, int nf), 
double slgg_plus(double x, int nf), 
double slp(double x,double y, int nf), 
double s2p(double x,double y, int nf), 
double slm(double x, double y, int nf), 
double s2m(double x,double y, int nf), 
double slp_plus(double x, int nf), 
double slm_plus(double x, int nf), 
double slgf_lim(double sp, double nf), 
double slfgJim(double sp, double nf), 
double s2ff_hm(double sp, double nf). 



24 



double s2fg_lim(double sp, double nf), 
double s2gf_lim(double sp, double nf), 
double s2gg_lim(double sp, double nf), 
double s2p Jim (double sp, double nf), 
double s2m_lim( double sp, double nf). 

Analytically computes or reads from the file the NLO weights for the evolution 
equation. These routines are grouped into 3 kinds: the sl,2xx routines calculate 
the regular weights, the sl,2xx_lim routines calculate the regular weights called 
at 1 and sl,2xx_plus do the weights that contain the plus-distributions. 

6. 5 alpha, c 

subroutines: 

double alpha(double tt, int nf), double alphae (double tt,int nf). 

Calculates LO, NLO and exact running coupling using corresponding pa- 
rameters from the input file. 

6. 6 init. c 
subroutines: 

double init q_uv( double xx), 
double initq_dv(double xx), 
double init_gl (double xx), 
double initq_ss(double xx), 
double init q_del( double xx), 
double initq_udbar (double xx). 



Sets initial values for all parton densities using the GRV98 input for LO and 
NLO densities from [^ . 

6.7 polylo.c 

subroutines: 

double Li2 (double x). 
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double Li3 (double x), 
double S12 (double x). 



Calculates these three polylogarithms using a fast routine with Bernouilli num- 
bers. . 

6.8 intpol.c 
subroutines: 

double int_q(int j, double xx,int it), 

double interpolate (double xx,double *xt, double *yt,int points). 

Interpolation routines used to calculate densities between grid points and for 
integration at the threshold. 

6. 9 evolver. c 

subroutines: 

evolver(int itl,int it2,int ic,int ib). 

The main routine that performs the evolution between thresholds for all den- 
sities. It updates the main density array q and the density derivatives array 
qp. 

6.10 thresh, c 

subroutines: 

int threshold(int what,int itt), 

int fdens4(double xx,int ittc, double *u, double *d, double *s), 

double light _charm( double xx,int ittc), 

double fcharm(double xx,int ittc), 

double fbottom(double xx,int ittc), 

double fsigma(doublc xx,int ittc), 

double fgluon(double xx,int ittc), 

double fcharm(double xx,int ittc), 

double fbottom (double xx,int ittc). 
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Threshold handhng routines to implement LO, NLO and NNLO matching con- 
ditions for light and heavy densities at the charm and bottom thresholds. The 
density routines are calls to convolution integrals that generate new densities 
for n/ + 1 flavors. 



6.11 a-coefs.c 
subroutines: 

double alqg(double z, double fs2,double hm2), 
double a2qq(double z, double fs2, double hm2), 
double a2qg(double z, double fs2, double hm2), 
double a2qqns(double z, double fs2, double hm2), 
double soft q( double z, double fs2, double hm2), 
double corq(double z, double fs2,double hm2), 
double a2gg(double z, double fs2, double hm2), 
double softg(double z, double fs2, double hm2), 
double corgi (double fs2, double hm2), 
double corg2(double z, double fs2,double hm2), 
double a2gq(double z, double fs2,double hm2). 



The OME routines used for NNLO threshold matching. These contain the 
formulae in Appendix B. 



6.12 loader, c 



subroutines: 

int loadOrd(int what). 

Functions to handle threshold datafile loading, saving and verification. This 
file allows one the ability to use previously computed density values at the 
threshold in a new computation. 
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6.13 quadrat, c 



subroutines: 

double qadr at (double *x, double a, double b, double (*fx) (double), double e[]), 
double lint(double *x, double (*fx) (double), double e[], double xO, double xn, 
double fO, double f2, double f3, double f5, double f6, double f7, double f9, 
double fl4, double hmin, double hmax, double re, double ae). 

Backup integration routine used as a check for the actual one used in the 
threshold integration. 



6.14 daind.c 
subroutines: 

double daind(double *x,double a,double b, double (*fun) (double), double eps,int 
key,int max). 

Main Gaussian integration routine, see [ P^ . 



6.15 integrands, c 
subroutines: 

inline double fcharm_integrand(double xl), 
inline double fgluonJntegrand(double xl), 
inline double fsigmaJntegrand(double xl), 
inline double us_integrand(double xl), 
inline double ds_integrand(double xl), 
inline double ss_integrand(double xl), 
inline double fbottom_integrand(double xl), 
inline double light _charmJntegrand( double xl). 



Functions containing integrands for the threshold integration. They use the 
density values and the coefficient functions from a-coefs.c to produce the in- 
tegrands that are then fed into the Gaussian integration program. 
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6.16 grids, c 



subroutines: 

int MakcXGrid(void), 

int MakcTGrid(void), 

int mcrgc(double *a, double *b,int na, int nb,char w), 
int check_grid(double *a,int n,char w), 
int MakeFortranGrid(int test_mode), 
double **allocate_real_matrix(int ur, int uc), 
void free_real_matrix(double **m,int ur). 



Subroutines for making (and also merging and verifying) the initial grids in 
X and and the final grids for Fortran-code compatible output. The grid 
merging is used to combine the evenly spaced grid generated automatically 
from the initial and final values with the premade grid containing several x 
and values for plotting and outputting the data. Two routines are added 
for deallocating memory. 

6.17 weights, c 

subroutines: 

int readWeights(int nf,int order), 
int dumpWeights(int nf,int order). 

Routines dealing with loading and saving computed NLO and NNLO weight 
tables to do a fast calculation on the same grids. LO weights are not saved as 
it is very fast to compute them every time. 

6.18 nnl-a-w.c 

subroutines: 

int anannlowgts(int nf,int loadWgts), 
int computeNNLOwgts(int nf), 
double nn_slff (double x,double y, int nf), 
double nn_s2ff(double x,double y, int nf), 
double nn_slfg(double x,double y, int nf), 
double nn_s2fg(double x,double y, int nf), 
double nn_slgf (double x,double y, int nf). 
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double nn_s2gf (double x,double y, int nf), 
double nn_slgg(double x, double y, int nf), 
double nn_s2gg(double x, double y, int nf), 
double nn_slff_plus(double x, int nf), 
double nn_slgg_plus(double x, int nf), 
double nn_slp(double x,double y, int nf), 
double nn_s2p(double x,double y, int nf), 
double nn_slm(double x, double y, int nf), 
double nn_s2m(double x, double y, int nf), 
double nn_slp_plus(double x, int nf), 
double nn_slm_plus(double x, int nf), 
double nn_slgf_lim(double sp, double nf), 
double nn_slfgJim(double sp, double nf), 
double nn_s2ff_lim(double sp, double nf), 
double nn_s2fg Jim (double sp, double nf), 
double nn_s2gfJim(double sp, double nf), 
double nn_s2ggJim(double sp, double nf), 
double nn_s2p_lim(double sp, double nf), 
double nn_s2m_lim(double sp, double nf). 

Analytically computes or reads from files the approximate NNLO weights for 
the evolution equations. Here the routines are grouped into three kinds: the 
nn_sl,2xx routines calculate the regular weights, the nn_sl,2xx_lim routines 
calculate the regular weights called at 1 and nn_sl,2xx_plus do the weights 
that contain the plus-distributions. 

6.19 wgplg.c 

subroutines: 

double wgplg(int n,int p, double x). 

The routines which calculate polylogarithms using the method from CERN- 
LIB [^. They are only used for the higher order polylogarithms because the 
routines for Li2, Li3 and S12 in polylo.c are faster. 
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7 Results 



The code can be used in several modes of operation. 

For all of them there is some optimum grid size in x and Q^. Internally, the 
grid with the sizes entered in the parameter file is merged with another grid 
(that is used for plotting the output data at the end), thus increasing the 
resulting grid size. This internal grid size contains all "popular" values, like 
X = 0.1, 0.01, 0.001 etc., and is 38 in and 64 in x. The corresponding values 
are located in file main.h (arrays xpr[] and q2pr[]). This grid is then merged 
with the automatically generated equidistant grid and the equal values are 
weeded out. Shown in the table are the resulting grid sizes as shown in the 
output file. The table uses the calculation for all fiavors as opposed to the 
GRV98-like (only three-fiavor) densities. In general, the evolution time grows 
quadratically in and linearly in nQ2. The numbers we give below arc for 
an alpha PC with a 21164 processor unit running at 500 MHz, 1 Gbyte of 
memory and rated at an Specfp = 20.4. 



order 


nx 




accuracy, digits 


time,sec 


LO 


162 


96 


5 


15 


NLO 


162 


96 


3 


113 


NNLO 


162 


96 


3 


385 


LO 


262 


136 


6 


31 


NLO 


262 


136 


5 


275 


NNLO 


262 


136 


5 


1021 


LO 


362 


136 


6 


44 


NLO 


362 


136 


6 


529 


NNLO 


362 


136 


5 


1537 



1. Set parameters to the following values to produce LO and NLO GRV98- 

style fixed three-fiavor densities for the whole range of (only parameters 
essential for this calculation are provided, the rest can be set to whatever one 
wishes since they control the form of the output and similar features, not the 
physically meaningful ones): 
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LO 




0.26 


Qinitial2 





Order 


1 


ThreeFlavorMode 


NLO 




0.40 


Qinitial2 


1 


Order 


1 


AlphaUseExact 


1 


ThreeFlavorMode 



2. To generate regular VFNS densities with all heavy flavors (both charm and 
bottom) one sets: 



LO 




0.26 


Qinitial2 





Order 





ThreeFlavorMode 


1 


AlphaDoSeparatcThreshold 





NNLOmultiOrderCHARM 


1 


DoBottomThreshold 


NLO 




0.40 


Qinitial2 


1 


Order 





ThreeFlavorMode 


1 


AlphaDoSeparateThreshold 


1 


AlphaUseExact 





NNLOmultiOrderCHARM 


1 


DoBottomThreshold 
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3. To generate VFNS densities involving proper order mixing at heavy thresh- 
olds with all heavy flavors (both charm and bottom) but without using NNLO 
weights (as done in our previous papers ||2^, [Q, ||3^) one sets: 



LO 




0.26 


Qinitial2 





Order 





ThreeFlavorMode 


1 


AlphaDoSeparateThreshold 


1 


DoBottomThreshold 


NLO 




0.40 


Qinitial2 


1 


Order 





ThreeFlavorMode 


1 


AlphaDoSeparateThreshold 


1 


AlphaUseExact 


1 


DoBottomThreshold 


NNLO 




0.40 


Qinitial2 


2 


Order 





ThreeFlavorMode 


1 


AlphaDoSeparateThreshold 


1 


AlphaUseExact 


1 


NNLOmultiOrderCHARM 


1 


DoBottomThreshold 


1 


NL04NNL0 



In this mode it is necessary to generate LO and NLO sets by running the 
program before running the NNLO set on the same grid! Those will be dumped 
in special data files 
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(agrv991o. BO. threshold, agrv991o.CH. threshold, 
agrv99nlo. BO. threshold, and agrv99nlo.CH. threshold) 
that will later be read for the NNLO calculation whenever 
NNLOmultiOrderCHARM=l . 

4. To generate VFNS densities involving proper order mixing at heavy thresh- 
olds with all heavy flavors (both charm and bottom) and using LO, NLO and 
NNLO (approximate) weights one sets: 



LO 




26 


Oinitial2 





Order 





ThreeFlavorMode 


1 


AlphaDoSeparateThreshold 


1 


DoBottomThreshold 


NLO 




n 4n 


Oinitial9 


1 

X 


OttI pt* 


n 


J- 111 C'C' J- Id V yji. ivi.wvJ.C' 




A Ifi n fi TjoSpt^;^ pT^ nrp^inoln 

li.llJllCX'l^yJkJy-'VJCXiL ObKjK^ J. Ill CollVjlLl 


1 


A 1 r>H Pi T TspF^jVji cf. 

-Til L/llC* \J Ovy J— J-A-C*^ U 


1 


DoBottomThreshold 


NNLO 




0.40 


Qinitial2 


2 


Order 





ThreeFlavorMode 


1 


AlphaDoSeparateThreshold 


1 


AlphaUseExact 


1 


NNLOmultiOrderCHARM 


1 


DoBottomThreshold 





NL04NNL0 



34 



In this mode it is also necessary to generate LO and NLO sets by running 
the program before nmning the NNLO set on the same grid! Those will be 
dumped in special data files 

(agrv99lo.BO.threshold, agrv99lo.CH.threshold, 
agrv99nlo. BO. threshold, and agrv99nlo.CH. threshold) 

that will later be read for NNLO calculation whenever 
NNLOmultiOrderCHARM= 1 . 

Program output is arranged in several forms. First, the default output in nor- 
mal readable form goes into resLO.dat, resNLO.dat or resNNLO.dat 
or for GRV98-mode into resL03.dat, resNL03.dat or resNNL03.dat 

depending upon the set calculation order. This file contains the input param- 
eters, calculation time and the columns of data versus and x for all densities 
(uv,dv,us,ds,ss,ch and bt, described in previous chapters). Here is the sample: 



========================== q2= 1960 ======================= 

Alpha(Q2= 1.96 GeV2)=0 . 318513 for nf=4 

x=0. 000010 

SI(x= 0.0000100) =3. 4695646e+00 GL(x= 0. 0000100) =1.3074834e+01 

UV(x= 0.0000100)=6.1120367e-03 DV(x= 0.0000100) =3. 7959190e-03 

US(x= 0.0000100=5. 9818110e-01 DS(x= 0.0000100) =5. 9988948e-01 
SS(x= 0. 0000100=5. 3175774e-01 

CH(x= 0.0000100)=0.0000000e+00 BT(x= 0. 0000100) =0.0000000e+00 

x=0. 000020 

SI(x= 0. 0000200) =3. 1153438e+00 GL(x= . 0000200) =1 . 1469641e+01 

UV(x= 0.0000200) =8. 2564168e-03 DV(x= 0. 0000200) =5. 1210226e-03 

US(x= 0. 0000200=5. 4149612e-01 DS(x= 0.0000200)=5.4372565e-01 
SS(x= 0.0000200)=4.6576141e-01 

CH(x= 0. 0000200) =0.0000000e+00 BT(x= 0. 0000200) =0.0000000e+00 



The above sample was produced with GraphAll=0 thus printing only values on 
a small grid with minimum = 1.96 GeV^ and not all values from minimum 
— 0.40 GeV^. For convenience, SI denotes singlet, GL gluon, UV and DV 
are valence densities u — u, d — u, US, DS, SS are of the q + q — 'E{nf)/nf kind 
and CH and BT are (c + c) /2 and {b + b) /2. 

For graphing purposes, the output also goes into several datafiles with names 
formed as g_densityORDER.dat where ORDER is LO, NLO or NNLO 
respectively e.g. g_glLO.dat or g_uvNNLO.dat. Those contains columns of 

the particular density with the first column being x or Q"^. depending upon 
GraphVsX parameter 0-Q^). Then the other parameter is varied across 

columns. Here is the piece of g_cpNLO.dat file. The first column contains 
the X value, the second is the charm density for = 1.96 GeV^ (where it is 
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zero) and then the charm density for = 2,3, .. GeV^: 



0.0000100000 
0.0000200000 
0.0000300000 
0.0000400000 
0.0000500000 
0.0000600000 
0.0000700000 
0.0000800000 
0.0000900000 



OOOOOOOOOOe+00 
OOOOOOOOOOe+00 
OOOOOOOOOOe+00 
OOOOOOOOOOe+00 
OOOOOOOOOOe+00 
OOOOOOOOOOe+00 
OOOOOOOOOOe+00 
OOOOOOOOOOe+00 
OOOOOOOOOOe+00 



1.0468420825e-02 
8.8559755303e-03 
8.0049032415e-03 
7.4395488912e-03 
7.0219632243e-03 
6.6940248888e-03 
6.4257535493e-03 
6.1998499386e-03 
6.0054517691e-03 



1.3022045486e-01 
1.0970003578e-01 
9.8909559249e-02 
9. 1758900075e-02 
8.6487022224e-02 
8.2352079221e-02 
7.8973989377e-02 
7.6132631291e-02 
7.3690103826e-02 



The above sample was produced with Graph VsX=l thus printing x, not 
values in the first column. The GraphAll=0 was also set, thus only nice values 
of X are used (0.00001, 0.00002, 0.00003, etc). 

Also, if the necessary option (DoFortran=l) is set the output also goes into 
the file suitable for reading by a GRV98-like Fortran program that interpolates 
the data points and makes parton density functions. This program is used in 
structure function calculations (the code is written in Fortran). The datafile 
format has eight columns with all densities on the fixed grid (hard-coded into 
the both evolution code and the interpolation program) in x and Q^. 

The sample follows: 

Information line: first 

+6.112E-03 +3.796E-03 +5.982E-01 +5.999E-01 +5.318E-01 +1.307E+01 

+6.128E-03 +3.806E-03 +6.084E-01 +6.101E-01 +5.419E-01 +1.333E+01 

+6.303E-03 +3.912E-03 +7.251E-01 +7.268E-01 +6.581E-01 +1.634E+01 

+6.440E-03 +3.996E-03 +8.260E-01 +8.278E-01 +7.586E-01 +1.901E+01 

+6.553E-03 +4.064E-03 +9.151E-01 +9.169E-01 +8.473E-01 +2.140E+01 

+6.649E-03 +4.122E-03 +9.949E-01 +9.967E-01 +9.269E-01 +2.357E+01 

+6.731E-03 +4.172E-03 +1.067E+00 +1.069E+00 +9.990E-01 +2.556E+01 

+6.804E-03 +4.216E-03 +1.134E+00 +1 . 135E+00 +1.065E+00 +2.740E+01 

+6.869E-03 +4.255E-03 +1.195E+00 +1 . 197E+00 +1.126E+00 +2.910E+01 

+6.927E-03 +4.291E-03 +1.252E+00 +1.253E+00 +1.183E+00 +3.070E+01 



Sample pictures of bottom densities are provided in and also below in 
Figs. 1 - 4. 
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8 Error code descriptions 

Program error code description: 



message 


filename 


refer to 


Threshold LO datafile 
is missing 


loader. c 


NNLO calculation with 
proper orders requires the 
datafile from a previous run 
in LO 


Threshold NLO 
datafile is missing 


loader.c 


NNLO calculation with 
proper orders requires the 
datafile from a previous run 
in NLO 


Wrong Multicharm 
factor 


main.c 


NNLOmultiOrderCHARM 
should be 1 or 


Wrong order factor 


several modules 


should be 0,1,2 for LO, 
NLO, NNLO 


File evolu- 
tion_parameters. input 
does not exist 


main.c 


find the file and put into 
working directory 


Wrong INI 
Q2: increase it! 


main.c 


order and initial are in- 
compatible 


Wrong INI 
Q2: decrease it! 


main.c 


order and initial are in- 
compatible 


Evolver: dont know 
how to proceed 


main.c 


wrong doBottom factor 


Wrong Alpha switch 
factor! 


main.c 


check AlphaDoSepa- 
rateThreshold value 


Wrong order factor 
while graphing 


main.c 


check Order to be 0,1,2 


Wrong graphing fac- 
tor 


main.c 


check GraphAU value to be 
0,1 


Wrong loadWgts fac- 
tor 


1-a-w.c, nl-a-w.c 


check loadWgts to be 0,1 
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9 Conclusions 



We have presented a multifunctional code for the direct x-space method of 
solving the spin-averaged evolution equations for parton densities. The dis- 
tinctive features of this code include analytic computation of the LO, NLO 
and NNLO weights, NNLO heavy flavor threshold matching and NNLO evo- 
lution. 

The code is very fast and accurate. For example for grid sizes not exceeding 
200 in and 150 in x the NLO calculation with full weighs computed for 
three values of and up to five decimal accuracy has a runtime well below 
200 seconds. Also it is the only code that does the proper NNLO evolution 
with NNLO heavy flavor matching conditions. 

The program is also easy to use and complete documentation is available. The 
code is well-tested both on specific test functions (e.g. see |T^) and on actual 
densities (e.g. see [^]) in all (LO, NLO and NNLO) orders. 
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A Appendix A 



Here we give tlie NNLO parametrizations of tlie splitting functions from p9 
Note that Lq = In^; and Li = ln(l — z). 

First the parametrizations for the non-singlet splitting functions Pj^s ^^6- 



P^2K(z)= 1185.229 (1 - z)-^ + 1365.458 5(1 - z) - 157.387 - 2741.42 z^ 

- 490.43 (1 - ^) + 67.00 + 10.005 Ll + 1.432 

+ Nf {-184.765 (1 - z)-^ - 184.289 5{1 - z) + 17.989 Lj + 355.636 

- 73.407 (l-2)Li + 11.491 + 1.928 Lj^} + PisU-^), 
pf^^^-(z)= 1174.348 (1 - z)~^ + 1286.799 5(1 - z) + 115.099 + 1581.05 Li 

+ 267.33 (1 - 2) - 127.65 Ll - 25.22 Lj^ + 1.432 
+ Nf {-183.718 (1 - z)-^ - 177.762 5(1 - z) + 11.999 + 397.546 z'^ 

+ 41.949 (1 - 2) - 1.477 Ll - 0.538 + Pns,2(^), (A.l) 



Pj^|+ (^)= 1183.762 (1 - z)-^ + 1347.032 5(1 - z) + 1047.590 Li - 843.884 z^ 

- 98.65 (1 - 2) - 33.71 Ll + 1.580 (L^ + ALl) 

+ A^j {-183.148 (1 - z)-^ - 174.402 5(1 - z) + 9.649 + 406.171 z^ 

+ 32.218 (1-2)+ 5.976 + 1.60 + Piiy^). 
pg+ (2)= 1182.774 (1 - 2);^ + 1351.088 5(1 - z) - 147.692 L\ - 2602.738 

- 170.11 + 148.47 Lo + 1.580 (L^ - 4L|j) 

+ iV/ {-183.931 (1 - z)-^ - 178.208 5(1 - z) - 89.941 Li + 218.482 z"^ 

+ 9.623 + 0.910 - 1.60 Lj^} + Pjiiy^) . (A.2) 
The parametrizations for P^^^'^ {z) and P^^{z) are 



pg^(2) = iVj {(1 - 2)(-1441.57 + 12603.59 z - 15450.01) + 7876.93 zLl 

- 4260.29 Lo - 229.27 Ll + 4.4075 L^} 
pg'|(2) = iVj {(1 - 2)(-704.67 z^ + 3310.32 + 2144.81 z - 244.68) 

+ 4490.81 z'^Lo + 42.875 Lq - 11.0165 Ll}, (A.3) 

and 



^psU('2)=^/ {(1 - ;z) (-229.497 Ll - 722.99 z^ + 2678.77- 560.20 z"^) 
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+ 2008.61 Lo + 998.15 - 3584/27 z'^Lo} + P^%(z), 
P^%B(z)^Nf {(1 - (73.845 Ll + 305.988 Li + 2063.19 z - 387.95 z'^) 
+ 1999.35 zLo - 732.68 Lq - 3584/27 z'^Lq} 
+ P^'U^), (A.4) 



with 



P§'l-2{z)^N'j {(1 - ^) (-7.282 Ll - 38.779 z^ + 32.022 ^ - 6.252 + 1.767 z''^) 
+ 7.453 L^} . (A.5) 

Next we show the parametrizations of the off-diagonal singlet splitting func- 
tions: 



PqtA{^)^Nf {-31.830 L\ + 1252.267 Li + 1999.89 z + 1722.47 + 1223.43 Ll 

- 1334.61^-^ - 896/3 z-^Lo} + Pq%{z), 
Pqg,B{^)^^f {19.428 L{ + 159.833 + 309.384 L\ + 2631.00 (1 - z) 



67.25 - 776.793 - 896/3 z-^Lo} + Pqg,2{z)^ (A.6) 



with 



Pgg'2{z)^Nf {-0.9085 Li - 35.803 Li - 128.023 + 200.929 (1 - z) 

+ 40.542 Lo + 3.284 z"^} , (A.7) 

and 

Pgf,Aiz)= 13.1212 Ll + 126.665 Lf + 308.536 Lj + 361.21 - 2113.45 Lq 

- 17.965 z-^Lo + Nf (2.4427 L^ + 27.763 Lf + 80.548 Lj 

- 227.135 - 151.04 L^ + 65.91 z'^Lo} + Pgf,2{z), 

pSA^)= -4-5108 L\ - 66.618 L\ - 231.535 L^ - 1224.22 (1 - ^) + 240.08 L^ 
+ 379.60 z~i(Lo + 4) + iVj (-1. 4028 Lf - 11.638 Lf + 164.963 Li 

- 1066.78 (1 - ^) - 182.08 L^ + 138.54 ^"^(Lq + 2)} 

+ Pl%{zl (A.8) 

with 



Pm,2{^)= {1-9361 L? + 11.178 Li + 11.632 - 15.145 (1 - ^) + 3.354 Lq 
- 2.133 z-^} . (A.9) 
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Last we show the parametrizations of the diagonal singlet splitting functions 



^isiAl-^)^ 2626.38 (1 - z)-'^ + 4424.168 S{1 - z) - 732.715 Lj - 20640.069 z 

- 15428.58 (1 - z^) - 15213.60 + 16700.88 z'^ + 2675.85 z'^Lq 
+ Nf {-415.71 (1 - z);^ - 548.569 S{1 - z) - 425.708 Li + 914.548 z'^ 

- 1122.86 - 444.21 + 376.98 z''^ + 157.18 z'^Lo} 

+ Pfkz), 

P^f j^{z)= 2678.22 (1 - z)-^ + 4590.570 5(1 - z) + 3748.934 Li + 60879.62 z 

- 35974.45 (1 + z^) + 2002.96 + 9762.09 z'^ + 2675.85 z'^Lo 
+ Nf {-412.00 (1 - z)-^ - 534.951 5(1 - z) + 62.630 + 801.90 

+ 1891.40 Lo + 813.78 + 1.360 z~^ + 157.18 z~^Lo} 

with 

PSU^)^^f {-16/9 (1 - z);^ + 6.4882 5(1 - z) + 37.6417 z^ - 72.926 z 
+ 32.349- 0.991 + 2.818^-^} . (A.ll) 
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B Appendix B 



Shown below are the renormahzed OME's used for threshold matching calcula- 
tions in NLO and NNLO (they correspond to the unrenormalized expressions 
given in Appendix C of and in Appendix A of 0]). All OME'S have been 
renormalized in the MS-scheme. 

In particular the renormalized coupling is presented in the above scheme 
for Uf + l light flavors. Here the heavy quark H = (c, b) is treated on the same 
footing as the light flavors and it is not decoupled from the running coupling 
in the VFNS approach. The (as/47r)^ coefficient in the heavy-quark OME 
is given by 



16 

'1 + z)\n.z 4 

"iz 



1 2 

111^^ + 



160 



9z 
+ (1 + ^) 

-Un^z 
3 



16 - 48z + 



448 



9 



1 + z) In^ z 

2 

m 
In — 

/i2 



+ A<dz + j In -2 



+ 1 



32 

'Yz 



32Si,2(l -z) + 161n2Li2(l - z) - 16C(2) Inz 



/32 32 A , 

- 8 + 82 + y J C(2) + [ 2 + IQz + y ^2 J 1^2 ^ 



56 88 



448 2 \ , 448 4 124 

z \mz 

/ 27^ 3 3 



z + 



1600 



(B.l) 



The as/47r and the (as/47r)^ coefficients of the heavy quark OME's yl^^ are 



is,(i)/m^ 



Tf 



-Aiz' + il 



nin^ 



(B.2) 



and 



78,(2) / 



-- |CfT/[(8 - 16z + 16^2) ln(l - ;2) 
-(4 - 8;^ + 16^2) In z - (2 - 8z)] 
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+CATf 



16 



-(8 - 16z + 16z^) ln(l -z)-{8 + 32z) In z 



— -A-32Z + 



124 



3z 
+{cFTf 



(8 - 16z + 16^^) [2 In z ln(l - - ln^(l - ^) + 2C(2)] 



-(4 - 8^ + 16^2) In^ z - 32^(1 - z) ln(l - ^) 
-(12 - 16z + 32^2) In z - 56 + 116^ - 80^^ 

+(7^7/ (16 + 32z + 32z2)[Li2(-^) +lnzln(l + z)] 

+ (8 - 16^ + 16z^) \n^{l - z) + {8 + 16z) In^ z 

( 352 \ 

+32<(2) + 32^(1 - z) ln(l - ^) - ( 8 + 64^ + —z^ \ In , 



160 

""97 



+ 16 - 200;^ + 



1744 



I VP? 



+CirT/{(l - 2z + 2z2)[8C(3) + - ln3(l - z) 
-81n(l - ^)Li2(l - ^) + 8C(2) Inz - 41n^ln2(l - z) 
+^ In^ ^ - 8 In zLi2(l - z) + 8Li3(l - ^) - 24Si,2(l - z)\ 

r 4 

^z^ -16C(2)ln;2 + -ln3;2+ 161n;2Li2(l + 32Si,2(l 

-(4 + 96z - 64^2)Li2(l - ^) - (4 - 48z + 40z^)C(2) 
-(8 + mz - 24^2) In ^ ln(l - + (4 + 8z - 12z^) W(l - 
-(1 + 12^ - 20^2) In^ z - {52z - 48^^) ln(l - z) 

-(16 + 18z + 48z^) In ^ + 26 - 82z + 80z^^ 

+CATfl^{l -2z + 2z2)[-l ln3(l - z) 

+8 ln(l - 2)Li2(l -z)- 8Lk{l - z)] + {1 + 2z + 2z^) 
x[-8C(2) ln(l + z) - 161n(l + 2)Li2(-z) -81nzln2(l + z) 
+41n2^1n(l + ;2) + 81n;2Li2(-;2) -8Li3(-;2) - 16Si,2(--2)] 

+(16 + 64z)[2Si,2(l - ^) + ln^Li2(l - z)\ - + \^ ^ 

+ (8 - 32^ + 16^2)C(3) - (16 + 64^)C(2) In^ + (16 + 16^^) 

(32 272 \ 

— + 12 + 64^- -^^MLi^^l 
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- f 12 + 48^ - —z^ + - )C(2) - 4^2 In ^ in(l - z) 
\ 2) 2>z J 

( 46 \ 

-(2 + 8z - lOz^) In^ (1 - ;2) + ( 2 + 8^ + y j In^ z 

/. 2M /-. N 172 1600 2\ , 

+(4 + 16z - 16^2) lii(l - - I y + -^-2 + 1 In -2 

448 4 628 6352 .l 

^ + zM, (B.3 

27^ 3 3 27 J' ^ ' 

respectively. Now we present the renormalized expressions for the heavy-quark 
loop contributions to the light-parton OME's denoted by A^i^h. The coeffi- 
cients of the {as/4:TT)'^ terms in Agg H and Agg^H are 



+ 



80 



9 \ 1 - ^ 



+ 



3 yi-z 

81 + ^2 



4 4 

---z + 25{l-z) 
8 88 



3 1 



+5(i-.)(^c(2)+?; 



\'n.z-\- —z 

z 9 9 

2 



1 + ^V2, 2 20, ^ 

H m 2; H mz 

1-^ I3 9 



+ -(1 - ^)ln^ + 



224 / 



1 



44 268 
21 \l-z) , ^27~^^ 



+ 



(B.4) 



and 



48,(2) /m- 



160 



16 16 8 

3; ~ y + 3^ 



+ 
+ 



9z 
4/2 
3 Vi 
1 /448 



160 128 ,32 32 16 , , , 



In 



m 



)8 / 10 \ 
ln2(l - + - f — - 10 + 8z j ln(l - z) 



21 \ z 



448 + 344^J 



(B.5) 



respectively. The coefficients of the Q;s/47r and [as/^T^Y terms in Agg^u are 
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Tf 



3 /i' 



(B.6) 



and 



48,(2) / 



+CATf 



8/ 1 ' 
3\ 



16 , ^ 16 2 

i(l + z)\nz-\ h4-4^ 

3z 3 



8 16 8 82 

\- -z z'^ 

32 3 3 3 



m 



/i2 



80 



16 



8(1 + z) In^ z + {2A + AOz) \nz + 64 - 32z 

16„ >, 80/ 1 \ 
_(l + ,),„, + _ _ 



^z^+A5{l-z 



+ 



184 232 152 184 



16 



Qz 9 



H z - — + —5il - z) 

9 9 3 ^ ' 



In 



m 



+CFTfl 3(1 + -2) In^ 2 + (6 + 10^) In^ 2 + (32 + 48^) In^ 



-- + 80 - 48z - 24^2 _ 155(1 - 

1(1 + -2) In^ z + ^(52 + 88z) \nz- ^z\n{l - z) 



1 

+ 27 



7 

2241 



9 
556 



628 + 548^ - 700^^ 



(B.7) 



respectively. 



The definitions for the polylogarithms Li„(z) and the Nielsen functions S„,p(2;), 
which appear in the above expressions, can be found in P7[] . 
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Figure Captions 



Fig. 1. The gluon density a;5'NNLo(4, a;, /i^) in tiie range 10 ^ < a; < 1 for 

IJ,^ = 2, 3, 4, 5, 10 and 20 in units of (GeV^)^, 
Fig. 2. The singlet density a;SNNLo(4, a;, /x^) in the range 10~^ < x < 1 for 

/i^ = 2, 3, 4, 5, 10 and 20 in units of (GeV^)^, 
Fig. 3. The nonsinglct quark density a;o"NNLo(4, a;, /i^)a, where a = {u + u)/2, 

in the range 10-5 < a; < 1 for /i^ = 2, 3, 4, 5, 10 and 20 in units of {GeV^, 
Fig. 4. The charm quark density a;cNNLo(4, x, ji"^) the range 10~^ < x < 1 for 

li^ = 1.96, 2, 3, 4, 5, 10 and 20 in units of (GeV^)^, 
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